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ABSTRACT 

Comprehensive numerical solutions of the steady state incompressible viscous flow over a 
three-dimensional backward-facing step up to Re = 800 are presented. The results are ob- 
tained by the least-squares finite element method (LSFEM) which is based on the velocity- 
pressure-vorticity formulation. The computed model is of the same size as that of Armaly’s 
experiment. Three dimensional phenomena are observed even at low Reynolds number. The 
calculated values of the primary reattachment length axe in good agreement with experi- 
mental results. 



1. INTRODUCTION 


It is well known that for backward-facing step flows the reattachment length obtained by 2D 
calculations cannot match with the experimental results for Re > 450. All the deviations 
of two-dimensional simulations from experimental results are mainly due to the growing 
effects of three-dimensionality as Reynolds number increases. In recent years, the devel- 
opment of numerical schemes in fluid dynamics has been concentrated on the solution of 
three-dimensional problems. Ku et al . 1 applied a pseudospectral matrix element method to 
simulate the same three-dimensional problem only up to Re = 450 by employing the primi- 
tive variable formulation. A simplified marker-and-cell finite-difference scheme is applied by 
Ikohagi and co-workers 2 for the same 3D problem in curvilinear coordinates. However, no 
comparisons have been made with either the experimental data or other simulation results, 
and thus no primary conclusions can be really made. 

In this paper, we report our results for the 3D backward-facing step flow, and compare 
with experimental results. We use the least-squares finite element method (LSFEM) 3,4,s . 
For incompressible viscous flows, the LSFEM is based on the first-order velocity-pressure- 
vorticity formulation. This method uses C° elements to discretize the equations and mini- 
mizes the £2 norm of the equation residuals, and results in a symmetric and positive- definite 
algebraic system which can be solved by simple matrix-free iterative methods. This method 
leads to a minimization problem, and the choice of elements are thus free of the limita- 
tion of the Ladyzhenskaya-Babuska-Brezzi(LBB) conditions 6 . The LSFEM is also free of 
any parameters. Furthermore, all unknown variables are solved in a fully coupled manner, 
and only simple physical boundary conditions are imposed, no artificial numerical boundary 
conditions need be devised. 

2. VELOCITY-PRESSURE- VORTICITY FORMULATION 

The mathematical formulation and numerical scheme of the LSFEM for 3D Navier-Stokes 
equations can be found in Jiang et al . 5 . For completeness, here we give a brief description of 
the method. The steady-state incompressible Navier-Stokes equations can be recast as the 
following ti on dim ensional first-order quasi-linear velocity-pressure- vorticity formulation 
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in which u denotes the velocity, p the pressure, f the body force, Re the Reynolds number, 
and D the flow domain. Here an independent vector, the vorticity tD = V x u, is introduced 
to yield the first-order form such that the C° element can be used. As proved by Jiang et 
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al. 5 , the following compatibility condition is required for three-dimensional cases to make 
the system elliptic and will not be detailed here, 


V • o> = 0. 


(4) 


The first-order elliptic system in Cartesian coordinates for three-dimensional problems can 
be expressed as 
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This first-order system has seven unknowns and eight equations. Jiang 5 proved that 
the system is determined by introducing a dummy variable ^ which satisfies the boundary 
condition <j> = 0 on the boundary. The gradient, V<f>, is then added to the equation of 
vorticity definition, and thus yields an elliptic and determined system with eight unknowns 
and eight equations. Taking divergence of this vorticity equation will lead to A<p = 0, thus 
(f> = 0 in f2. Note that the introduction of a dummy variable <f) is purely for the purpose of 
proving the ellipticity of the system (5). In real calculation, no extra unknown is needed. 

As shown in Jiang et al. 5 , the 3D incompressible Navier-Stokes equations must have 
three boundary conditions on each boundary. For the backward-facing step flow problem, 
the following boundary conditions are considered : 

(a) u = t) = n; = 0on the wall; 

(b) u = specified , v — w = 0 at the well-developed inflow; 

(c) v = w = 0, p = constant (reference) at the outflow; 

(d) u> x = u) y = 0, w = 0 at the symmetric plane. 
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Note that no derivatives are involved in the boundary conditions, only the simple physical 
boundary conditions are imposed, and vorticity boundary conditions are used only at the 
symmetric plane. 

3. NUMERICAL SCHEME 


First, the Newton’s scheme is employed to linearize the quasi-linear system, thus we have, 
for example, 


u 


du 

dx 


u 


du 

dx 


+ <rJ 




where ( )° denotes the evaluation of the variable from previous iteration level, 
linearized first-order equations can be written as a standard first-order system 


( 6 ) 


Then the 


L u = f in 


( 7 ) 


where L is a first-order partial differential operator : 

. du du du 

L u = Ai — + A 2 - 5 - + A$ — + Au ( 8 ) 

dx dy dz 

The coefficient matrices in the general form of the first-order system for three-dimensional 
problems are 
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and the force and state variable vectors are 


0 

“”(£)” + *° (ft)” + «’"(£)” 
u°(U)° + v°(U)° + w °& 

+»°(%r 

0 
0 
0 
0 

The linearized first-order equations are solved by the LSFEM 3 which results in a sym- 
metric and positive-definite algebraic equation. A Jacobi preconditioned conjugate gradient 
method is employed to solve the linear algebraic equations. In the conjugate gradient method, 
the major computation is the multiplication of the global matrix with the global vector, and 
this can be done in an element-by-element manner without forming the global matrix 7 . To 
further save the storage, the current algorithm does not even form the element matrices. 
We directly calculate the product of the element matrix and the element vector, and the 
J acobi preconditioner can also be easily formed at the same time. In this way, we store only 
several global vectors, and the derivatives of the shape functions at Gaussian points for each 
element. 

4. NUMERICAL RESULTS 

The geometry and boundary conditions for the 3D model are shown in Figure 1. Due to 
the symmetry, only half of the domain is considered. The step has a height S=4.9 mm, and 
the inlet boundary is located 3.5 step heights upstream of this step. The channel downstream 
of the step has a height of 10.1 mm that provides an expansion ratio of 1 : 1.9423, and the 
half-span W = 90 mm that provides an aspect ratio for the test channel is 18 : 1.01. The 
length L measured from the step to the end of the calculation domain is 45 step heights, which 
is 3.11 times the maximum experimentally measured reattachment length of the primary 
recirculation region for the Renolds numbers interested. The Reynolds number Re = UDju 
is based on the hydraulic diameter (D = 10.4mm) of the inlet channel, and the average 
velocity is two-thirds of the maximum inlet velocity (normalized to unity) on the mid-span 
plane. The various Reynolds numbers are obtained by varying the kinematic viscosity v. 
The velocity profiles obtained by solving the 2D Poisson’s equation for a well- developed flow 
are imposed as the inflow boundary conditions. The outflow boundary conditions are defined 
such that a parallel flow and a constant pressure field exist. 

The computation was performed on the mesh with 54400 nonuniform trilinear elements 
(6 x 16 x 20 for inlet channel and 82 x 32 x 20 for the test channel). The smallest element 
has the size of 0.1251, 0.05714 and 0.421 steps in x,y and z direction individually, while the 
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largest element has the size of 3.0, 0.06938 and 1.68367 steps respectively. Figure 2 shows 
the nonuniform mesh which has more elements close to the sidewall, the floor and the roof 
in the test channel (see Figure 2). The solutions of the Stokes problem are taken as initial 
guesses for the case of Re = 100, and the ’’converged solutions are then used as initial 
guesses for higher Reynolds numbers. The ’’converged” solutions are based on those whose 
I 2 norm of the residuals less than 1.0 x 10 -4 . The problem is solved using about 5 M words 
of the memory on a CRAY-YMP. Simulation of the three-dimensional model is performed 
from Re = 100 to 800. 

S panwise flow structure 

The velocity profiles at three spanwise locations for Re = 277 and 800 are shown in 
Figures 3 and 5. At Re - 277, there is a slight change in the spanwise flow structure and 
reattachment length. At Re - 800, the velocity profiles change significantly in spanwise 
direction, and the reattachment length increases rapidly as moving toward the mid-plane. 
Figures 4 and 6 depict the pressure contours at these two Reynolds numbers. Again, the 
change in spanwise pressure distribution is more obvious as Reynolds number increases. Most 
researchers predicted the reattachment length well up to Re < 450 by their 2D simulations, 
and thought that the two-dimensional phenomena were maintained until Re « 450. Numer- 
ical results from the present method, however, show that the three-dimensionality are quite 
significant even at low Reynolds number (e.g. Re= 277). 

Figure 7 illustrates the contours of vorticity u> x from z = 6.5 cm to the sidewall ( z = 
9.0 cm) at Re = 277. The dashed lines indicate the negative contour values. As shown in the 
Figures, the vortex is stronger near the sidewall and has negligible influence on the region near 
the mid- plane. The study of velocity vectors along the span at different downstream locations 
provides better view of the three-dimensional phenomena, see Figure 8. The behaviour of the 
inward flow toward the mid-span at the top roof and the outward flow toward the sidewall 
(zjW = 1.0) at the channel floor contributes the three-dimensional phenomena. This three- 
dimensionality aro un d the comer of the sidewall and the step can also be found at lower 
Reynolds numbers but less significant. 

A series of plots of cross-flow velocity vectors and contours of stream wise vorticity, u> x , at 
different downstream locations are shown in Figures 9 - 14. At Re — 389, there is a counter- 
clockwise (negative) vortex on the comer of sidewall and channel floor. The cross-flow 
velocity vectors shown in Figure 10 depict that there is a weaker clockwise (positive) vortex 
close to the ch ann el ceiling as flow moving downstream. The size of this positive vortex grows 
in the streamwise direction. Both vortexes become stronger as Reynolds number increases, 
and persist further downstream with decreasing strength. 

It is found that except in the inlet channel and in the region which is very fax downstream 
of the step, the three-dimensionality is significant at the downstream, in the vicinity of the 
step. A study of spanwise flow structure provides further details. The spanwise distribution 
of velocity profiles for streamwise velocity, u, at various x-locations and at a fixed t/- position 
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for Re = 648 axe shown in Figure 15. At y = 7.5 mm, all velocity profiles close to the 
mid-plane axe basically two-dimensional. The negative velocities in the figure indicate that 
a second separation bubble occurs on the ceiling, and its thickness grows as close to the 
sidewall. At y = 2.35 mm, the x-component velocity, u, first increases rapidly in about the 
thickness of boundary layer, then drops and resumes two-dimensional flow from z/W = 0.5 
to the mid-plane ( zjW — 0.0). 

Figure 16 demonstrates the three-dimensionality by depicting the spanwise distribution 
of the reattachment length for the primary sepaxated-flow region. The numerical data of 
spanwise reattachment length vs. various Reynolds numbers axe listed in Table 1. The 
reattachment length is pretty much constant close to the mid-plane and decreases as moving 
toward the sidewall. It is interesting to note that next to the sidewall, the reattachment 
length increases rapidly. This phenomenon might be due to the interaction of primary 
recirculation vortex and the comer vortex between the sidewall and floor as shown in, for 
example, Figures 9 and 10. Figure 17 shows the computed reattachment length of the 
primary recirculation zone, and compares with Armaly et a/.’s 8 experimental results. The 
corresponding pointwise data obtained by the LSFEM and the experiment axe given in Table 
2. Since no tabular results are given, here the cited data from experiment axe obtained by 
optically digitizing Figure 14 in Armaly et a/.’s paper. 

Table 2 shows that in the Re range for which most simulations fail to predict the reattach- 
ment length because of the three-dimensional phenomenon (Re = 450 ~ 800), the calculated 
results by the LSFEM agree very well with the experimental results up to Re = 800. 

As Reynolds number increases, an additional sepaxated-flow region occurs near the chan- 
nel ceiling. Figure 18 illustrates the spanwise detachment (x^/S) and reattachment (x 5 /S) 
lines of this second eddy. The present results show that as Reynolds number increases, 
this upper-wall eddy propagates toward the mid-plane with its length decreasing toward the 
mid-plane. For example, as Reynolds number increases from 600 to 800, its length changes 
from 17.5 to 22 step heights at z/W = 0.97708, and changes from 0.2 to 9.36 step heights 
at z/W = 0.76. The experimentally observed upper-wall eddy at the mid-plane was not 
observed in the current simulation. The spanwise variation of detachment and reattachment 
length of the second eddy for various Reynolds numbers are given in Tables 3 and 4. 

5. CONCLUSIONS 

The steady-state three-dimensional backward-facing step problem is simulated using the 
least-squares finite element method. The computed spanwise flow structure cleaxly depicts 
the three-dimensionality. The prediction of primary reattachment length axe in good agree- 
ment with experimental results. Further developments are under way for solving time- 
dependent problems. 
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Table 1. Span wise Distribution of Primary Reattachment Length (si/5) 


z/W 

i2e=389 

.Re=500 

O 

O 

CO 

II 

ft} 

i2e=648 

O 

o 

II 

CD 

O 

O 

00 

II 

0.97708 

4.7873 

4.7598 

4.8502 

4.7166 

4.8231 

4.7052 

0.95243 

4.5224 

4.5636 

4.5934 

4.5244 

4.5810 

4.5098 

0.92591 

4.4467 

4.4700 

4.4290 

4.3849 

4.3936 

4.3248 

0.89739 

4.8376 

4.7809 

4.6438 

4.5963 

4.5364 

4.4150 

0.86670 

5.6686 

5.5985 

5.4196 

5.3551 

5.2089 

5.0171 

0.83370 

6.5050 

6.6123 

6.5767 

6.5552 

6.4553 

6.2616 

0.79819 

6.9495 

7.2317 

7.3644 

7.3418 

7.4005 

7.4756 

0.76000 

7.2835 

7.6769 

7.8567 

7.7813 

7.8771 

7.9793 

0.71892 

7.6081 

8.2233 

8.4387 

8.2980 

8.4029 

8.4371 

0.67472 

7.8156 

8.8039 

9.1300 

8.9540 

9.0833 

8.9984 

0.62719 

7.9101 

9.3172 

9.8642 

9.7360 

9.8871 

9.6695 

0.57605 

7.9336 

9.6589 

10.5745 

10.6628 

10.8592 

10.6415 

0.52104 

7.9291 

9.7898 

11.0796 

11.5275 

11.8613 

11.9200 

0.46187 

7.9188 

9.8100 

11.3667 

12.1405 

12.7017 

13.3638 

0.39823 

7.9056 

9.7822 

11.4917 

12.4468 

13.2472 

14.7292 

0.32976 

7.8933 

9.7539 

11.5069 

12.4740 

13.4493 

15.5563 

0.25611 

7.8839 

9.7385 

11.4735 

12.3895 

13.4234 

15.6984 

0.17689 

7.8780 

9.7272 

11.4279 

12.3410 

13.3384 

15.4442 

0.09167 

7.8748 

9.7209 

11.4000 

12.3408 

13.2772 

15.1296 

0.00000 

7.8740 

9.7186 

11.3926 

12.3450 

13.2504 

15.0204 


Table 2. Primary Reattachment Length (xi/S) 


Re 

389 

500 

600 

648 

700 

800 

LSFEM 

Exp. 

7.874 

8.67 

9.719 

10.21 

11.393 

11.40 

12.345 

12.36 

13.253 

13.10 

15.020 

14.45 
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Table 3. Detachment Length (x 4 /S) of Secondary Eddy at Upper Wall 


z/W 

i?e=389 

Re=500 

7Ze=600 

Re=648 

Re=700 

O 

o 

oo 

II 

<0 

ft? 

0.97708 

3.0831 

3.0024 

2.8182 

2.8258 

2.7367 

2.6701 

0.95243 

4.2783 

4.1200 

3.8681 

3.8440 

3.7324 

3.6367 

0.92591 

5.6721 

5.2026 

4.9072 

4.7961 

4.6904 

4.5429 

0.89739 

7.8333 

6.1612 

5.7713 

5.5829 

5.5169 

5.3662 

0.86670 

— 

7.3190 

6.6515 

6.3569 

6.3083 

6.1724 

0.83370 

— 

- 

7.6409 

7.2013 

7.1027 

6.9149 

0.79819 

— 

- 

9.0375 

8.2846 

8.0603 

7.6900 

0.76000 

— 

- 

11.3750 

9.6379 

9.3417 

8.7960 

0.71892 

- 

— 

- 

10.7500 

10.5547 

10.0307 

0.67472 

- 


- 

- 

11.9552 

11.0541 

0.62719 

— 

- 

- 

- 

- 

12.2672 

0.57605 

“ 

- 

- 

- 

- 

14.1967 


Table 4. Reatachment Length (x 5 /S) of Secondary Eddy at Upper Wall 


z/W 

0.97708 

0.95243 

0.92591 

0.89739 

0.86670 

0.83370 

0.79819 

0.76000 

0.71892 

0.67472 

0.62719 

0.57605 


Re = 389 
13.7775 
13.3167 
12.2407 
9.3224 


Re=500 

17.2220 

17.0891 

16.7674 

15.3696 

12.8966 


J?e=600 

20.3110 

20.0150 

19.6285 

18.3668 

16.5111 

14.7052 

13.3818 

11.5714 


i2e=648 

21.1037 

20.9073 

20.9271 

19.9921 

17.8685 

15.9967 

14.9380 

14.5169 

14.1530 


i£e=700 

22.7215 

22.5413 

22.4388 

21.4839 

19.5302 

17.7279 

16.7597 

16.1456 

15.7025 

14.8041 


jRe=800 

24.6619 

24.2501 

24.6982 

24.0531 

22.0990 

19.8740 

18.4722 

18.1593 

18.3837 

18.5512 

18.1801 

16.7677 
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Figure 2. Nonuniform mesh ( xy plane and xz plane ). 
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x/S=0 

















(c) z/W=0.0 (mid-plane) 
Figure 5. Velocity profiles for Re=800. 
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Figure 6. Pressure contours for Re=800. 













Figure 7. COx contours for Re=277 at x/S= (a) 6.25, 
(b) 9.0, (c) 12.25. 


Figure 8. Velocity vectors for Re=277 at x/S= (a) 6.25, 
(b) 9.0, (c) 12.25. 
















Figure 1 1. cox contours for Re=648 at x/S= (a) 6.25, Figure 12. Velocity vectors for Re=648 at x/S= (a) 6.25, 
(b) 9.0, (c) 12.25. 0>) 9 0 > ( c ) 12 25 - 
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Figure 17. Comparison of experimental and simulation results for 
primary reattachment length up to Re=800. 
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